Lecture 10: Spatial Data in R

class: center, middle, inverse, title-slide

.title[
# Lecture 10: Spatial Data in R
]
.subtitle[
## Vector Data Manipulation
]
.author[
### James Sears <br> AFRE 891 SS 24 <br> Michigan State University
]
.date[
### .small[<br> *Parts of these slides are adapted from <a href="https://github.com/msu-econ-data-analytics/course-materials">“Advanced Data Analytics”</a> by Nick Hagerty and <a href="https://github.com/dlab-berkeley/R-Geospatial-Fundamentals">“R Geospatial Fundamentals”</a> by the UC Berkeley D-Lab used under <a href="https://creativecommons.org/licenses/by-nc-sa/4.0">CC BY-NC-SA 4.0</a>.]
]

---

# Table of Contents: Today

**Part 2: Vector Data Manipulation**

1. [Spatial Queries: Measurement](#measure)

1. [Spatial Queries: Relationship](#relate)

1. [Geometric Operations](#geom)

1. [Spatial Joins](#join)

---

# Table of Contents: Next Time

**Part 3: Raster Data**

1. Common Raster Data Sources

1. Raster Operations

1. Combining Raster and Vector Data

---
# Prologue

To start, load in some packages/data from Part 1:

```r
if (!require("pacman")) install.packages("pacman")
pacman::p_load(sf, tidyverse, tmap, units)

counties<- st_read("data/MichiganCounties/MichiganCounties.shp")
```

```
## Reading layer `MichiganCounties' from data source 
##   `F:\OneDrive - Michigan State University\Teaching\MSU 2023-2024\AFRE 891 SS24\Lecture-Slides\10-Spatial\data\MichiganCounties\MichiganCounties.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.41829 ymin: 41.69613 xmax: -82.41348 ymax: 48.26269
## Geodetic CRS:  WGS 84
```

```r
ingham <- filter(counties, NAME == "Ingham")
```
---
# Preamble

Let's also define a custom ggplot map theme

```r
# add ggplot map theme
maptheme <- theme(
  panel.background = element_rect(fill = NA),
  # panel.border = element_rect(fill = NA, color = "grey75"),
  axis.text.x=element_blank(), #remove x axis labels
  axis.ticks.x=element_blank(), #remove x axis ticks
  axis.text.y=element_blank(),  #remove y axis labels
  axis.ticks.y=element_blank(),  #remove y axis ticks
  panel.grid.major = element_blank(),
  panel.grid.minor = element_blank(),
  legend.key = element_blank())
```

---
class: inverse, middle
name: measure

# Spatial Queries: Measurement

---

# Spatial Queries

.hi-blue[Spatial queries] ask questions of our data and return data metrics, subsets, or new data objects.

The two basic types of spatial queries are

.hi-medgrn[1\. Spatial Measurement]
- Return .hi-medgrn[continuous, numerical] answers

- .hi-slate[Area:] How many square meters is Lake Superior?

- .hi-slate[Length:] What is the length of I-96 in Michigan?

- .hi-slate[Distance:] How far is East Lansing from Traverse City?

---

# Spatial Queries

.hi-blue[Spatial queries] ask questions of our data and return data metrics, subsets, or new data objects.

.hi-purple[2\. Spatial Relationships]
- Return .hi-purple[TRUE or FALSE] (binary predicates)

- .hi-slate[Intersects, Overlaps, Touches, Disjoint, Crosses]
  - .hi-slate[Contains, Covered by, Covers, Within]
  - .hi-slate[Equal]
  
--

- Or the .hi-purple[set of matching features]

- .hi-slate[Intersection]
  - .hi-slate[Difference]
  - .hi-slate[Union]
  - .hi-slate[Crop]

---

# Spatial Queries

Note: sf can do a [*lot*](https://osf.io/an6b5/download) of additional spatial manipulations

---
# Spatial Queries

Note: sf can do a [*lot*](https://osf.io/an6b5/download) of additional spatial manipulations

<img src="data:image/png;base64,#images/sf-2.png" width = "88%"/>
---
# Calculating Area

```r
st_area(ingham)
```

```
## 1449829600 [m^2]
```

Too big!

```r
st_area(ingham)/1000000
```

```
## 1449.83 [m^2]
```

Wait, the units are wrong!

```r
units::set_units(st_area(ingham), km^2)
```

```
## 1449.83 [km^2]
```

[Resource: Measurement units in R](https://cran.r-project.org/web/packages/units/vignettes/measurement_units_in_R.html)
---
# Calculating Area

Next, let's add an `area` variable to our dataframe for each county:

```r
counties2 <- mutate(counties,
                   area = units::set_units(st_area(counties), km^2))
mean(counties2$area)
```

```
## 1818.946 [km^2]
```
--

```r
head(counties2$area)
```

```
## Units: [km^2]
## [1] 1795.1265 2420.5733 2178.4515 1536.1768 1356.7017  951.2289
```
---
# Calculating Area: Impact of CRS

Recall that spatial measurements can .hi-mdgrn[differ greatly] depending on the .hi-medgrn[CRS/projection we use!]

```r
# Calculate area using data in WGS84 CRS (4326)
counties2$area_wgs84 <- st_transform(counties, 4326) %>% st_area() %>% set_units(km^2)
  
# Calculate area using data in Web Mercator CRS (3857)
counties2$area_web <- st_transform(counties, 3857) %>% st_area() %>% set_units(km^2)
  
# Calculate are using data in UTM Zone 16N, NAD83 (26916)
counties2$area_utm16 <- st_transform(counties, 26916) %>% st_area() %>% set_units(km^2)
  
# Take a look at the results
counties2 %>% select(starts_with("area")) %>% st_drop_geometry() %>% head()
```

```
##               area       area_wgs84        area_web       area_utm16
## 1 1795.1265 [km^2] 1795.1265 [km^2] 3559.229 [km^2] 1800.7854 [km^2]
## 2 2420.5733 [km^2] 2420.5733 [km^2] 5102.846 [km^2] 2424.9251 [km^2]
## 3 2178.4515 [km^2] 2178.4515 [km^2] 4028.376 [km^2] 2180.7878 [km^2]
## 4 1536.1768 [km^2] 1536.1768 [km^2] 3082.999 [km^2] 1541.0568 [km^2]
## 5 1356.7017 [km^2] 1356.7017 [km^2] 2719.403 [km^2] 1359.3591 [km^2]
## 6  951.2289 [km^2]  951.2289 [km^2] 1846.495 [km^2]  953.8298 [km^2]
```

---
# Calculating Area: Impact of CRS

.hi-slate[Your choice of CRS is absolutely critical to accurate calculations.]

* *.hi-blue[WGS84] is a geographic (unprojected) CRS in degrees... but `sf` calculates area using spherical geometry!

* .hi-medgrn[UTM16N] is optimized for most of MI. Calculations outside the zone will be increasingly distorted as we move farther away from the zone
* .hi-purple[Web Mercator] preserves shape and explicitly .hi-purple[distorts area, don't use it for area calculations!]

---
# Calculating Length
`st_length()` calculates .hi-medgrn[feature lengths] in a spatial dataframe.

Let's first load in some data on railroad locations in Michigan.

.font80[

```r
#https://gis-michigan.opendata.arcgis.com/datasets/Michigan::railroads-v17a/explore?location=42.470352%2C-84.644141%2C7.94
  ## load in the railroads data
rail <- st_read("data/Railroads/")
```

```
## Reading layer `Railroads_(v17a)' from data source 
##   `F:\OneDrive - Michigan State University\Teaching\MSU 2023-2024\AFRE 891 SS24\Lecture-Slides\10-Spatial\data\Railroads' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 650 features and 10 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -90.21952 ymin: 41.7071 xmax: -82.42273 ymax: 46.79023
## Geodetic CRS:  WGS 84
```
]
---
# Calculating Length

Let's first load in some data on railroad locations in Michigan.

.font80[

```r
st_crs(rail)
```

```
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
```
]

.hi-slate[Q: should we reproject the railroad data before continuing?]

---
# Calculating Length
`st_length()` calculates .hi-medgrn[feature lengths] in a spatial dataframe.

```r
rail <- mutate(rail,
               len = st_length(rail)) 
head(rail$len)
```

```
## Units: [m]
## [1]   120.1038   142.0916 10207.5456 11527.4791 26872.0126   455.0065
```

Note the output units! (meter)

This mostly matches the (unitless) `LENGTH` variable that was already in the dataframe

---
# Calculating Length
Use `set_units()` from the .hi-slate[units] package to change to preferred units:

```r
rail <- mutate(rail,
               len_mi = units::set_units(st_length(rail), mi),
               len_km = units::set_units(st_length(rail), km)
               )

select(rail, starts_with("len")) %>% st_drop_geometry() %>% head()
```

```
##       LENGTH            len           len_mi          len_km
## 1   120.3320   120.1038 [m]  0.07462907 [mi]  0.1201038 [km]
## 2   142.1354   142.0916 [m]  0.08829164 [mi]  0.1420916 [km]
## 3 10220.1411 10207.5456 [m]  6.34267476 [mi] 10.2075456 [km]
## 4 11551.5451 11527.4791 [m]  7.16284343 [mi] 11.5274791 [km]
## 5 26889.6901 26872.0126 [m] 16.69749452 [mi] 26.8720126 [km]
## 6   456.0733   455.0065 [m]  0.28272792 [mi]  0.4550065 [km]
```

---
#Calculating Distance

How far are the railroads from the centroid of MSU's campus?

Let's create an sf object from point coordinates we [looked up on latlong.net](https://www.latlong.net/place/michigan-state-university-east-lansing-mi-usa-4710.html)

```r
msu <- data.frame(latitude = 42.701847,
                  longitude = -84.482170) %>% # build dataframe
  st_as_sf(coords = c("longitude", "latitude"), # specify variables with coordinates
           crs = 4326) # set WGS84 CRS
```
---
#Calculating Distance

Plotting on top of the Ingham County polygon:

```r
ggplot() +
  geom_sf(data = filter(counties, NAME == "Ingham")) +
  geom_sf(data = msu) + 
  maptheme
```

---
#Calculating Distance

`st_distance(X, Y)` calculates .hi-blue[spherical distances] by default

```r
rail_msu_dist <- st_distance(rail, msu)
```

---
#Calculating Distance

What does `st_distance()` return?

```r
class(rail_msu_dist)
```

```
## [1] "units"
```

```r
dim(rail_msu_dist)
```

```
## [1] 650   1
```

```r
head(rail_msu_dist)
```

```
## Units: [m]
##          [,1]
## [1,] 577183.6
## [2,] 548383.0
## [3,] 548383.0
## [4,] 577183.6
## [5,] 523377.8
## [6,] 487639.0
```

---
#Calculating Distance

By default, `st_distance()` returns a `unit` object with

* One row per feature in X (650 rail lines)

* One column per feature in Y (1 point for MSU)

---
#Calculating Distance
What about the distance to the .hi-medgrn[closest rail line] for our `counties` object?

```r
rail_cnty_dist <- st_distance(counties, rail) %>% # unit object, 1 row per county
  as.data.frame() %>% # convert to data frame
  rowwise() %>% # group by row 
  summarise(min = min(c_across(everything()))) # get row min across all columns

counties$rail_dist <- rail_cnty_dist$min # add as column to counties 
```

Which counties don't have a rail line in them?

```r
counties$NAME[which(counties$rail_dist > 0)]
```

```
## Error in Ops.units(counties$rail_dist, 0): both operands of the expression should be "units" objects
```
That doesn't work since our variable has units!

---
#Calculating Distance

.hi-medgrn[Solution:] either need to compare to a units object or change the variable to numeric.

```r
counties$NAME[which(counties$rail_dist > set_units(0, m))]
```

```
##  [1] "Benzie"       "Cheboygan"    "Gladwin"      "Keweenaw"     "Leelanau"    
##  [6] "Mecosta"      "Montcalm"     "Montmorency"  "Oceana"       "Oscoda"      
## [11] "Presque Isle"
```

---
#Calculating Distance

```r
ggplot() +
  geom_sf(aes(fill = as.numeric(rail_dist)), data = counties) +
  geom_sf(data = rail, color = "forestgreen") + 
  maptheme
```

---
class: inverse, middle
name: relate

# Spatial Queries: Relationship

---
# Spatial Queries: Relationship

Often we want to understand how spatial objects relate to each other.

.hi-purple[Binary predicates] tell us whether a specified relationship holds between geometries, returning a matrix of .hi-purple[TRUE or FALSE]

Common binary predicates include
* `st_intersects()` (most general)
* `st_within()`
* `st_contains()` (inverse of `st_within()`)

We can also get more specific about the type of geometric confirmation
* `st_overlaps()`, `st_touches()`, `st_disjoint()`, `st_crosses()`
* `st_covered_by()`, `st_covers()`
* `st_equal()`

---
# Spatial Queries: Relationship

Let's load in data on [protected wellhead areas](https://gis-michigan.opendata.arcgis.com/datasets/egle::wellhead-protection-areas-whpa/explore?location=42.680891%2C-84.545043%2C10.84) in MI to investigate

```r
whpa <- st_read("data/Wellhead_Protection_Areas_(WHPA)") %>%
  # fix invalid polygons
  st_make_valid()
```

```
## Reading layer `Wellhead_Protection_Areas_(WHPA)' from data source 
##   `F:\OneDrive - Michigan State University\Teaching\MSU 2023-2024\AFRE 891 SS24\Lecture-Slides\10-Spatial\data\Wellhead_Protection_Areas_(WHPA)' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2837 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -90.22974 ymin: 41.71514 xmax: -82.52524 ymax: 47.46764
## Geodetic CRS:  WGS 84
```

```r
# Filter county shapefile to just Ingham
ingham <- filter(counties, NAME == "Ingham")
```
---
# Spatial Queries: Relationship

And plot Ingham county with the wellhead protection areas...

... but first verify that they're on the same CRS!

```r
# Verify that CRS is the same as the rail line
stopifnot(st_crs(ingham) == st_crs(whpa))
```

---
# Spatial Queries: Relationship

And *actually* plotting Ingham county with the wellhead protection areas

```r
tm_shape(whpa) + tm_polygons() +
  tm_shape(ingham) + tm_polygons(col = "green", alpha = 0.5)
```

---
# Spatial Queries: Relationship
Some protected areas are completely .hi-blue[within] the county, some .hi-green[overlap], and some are completely .hi-purple[outside].

Our choice of geometric predicate will therefore influence *which* wellhead areas we select

Let's work through them and see how they differ.

---

# Relationship Predicates

`st_intersects()` is the .hi-medgrn[most general] of confirmation functions

* returns `TRUE` if x and y geometry share .hi-medgrn[any space],  `FALSE` if .hi-medgrn[absolutely no shared space]

We could also be more specific with the form of space that's shared

* `st_crosses()` looks for .hi-purple[any commonalities] between x and y.

* `st_touches()` looks if x and y .hi-green[share a common point] but .hi-green[interiors don't intersect]

* `st_overlaps()` looks if x and y .hi-pink[share space], are .hi-pink[of the same dimension], but .hi-pink[aren't completely contained by one another]

---

# Relationship Predicates

We could also be more specific with the form of space that's shared

* `st_within()` looks if x is .hi-blue[completely within] y

* `st_is_within_distance()` looks if x is .hi-red[within a specified distance] of y

* `st_contains()` looks if x .hi-dkorange[completely contains] y.

* `st_covered_by()` looks if x is .hi-slate[completely covered by] y.

* `st_covers()` looks if x .hi-yellow[completely covers] y.

* `st_disjoint()` looks if x .hi-grey[doesn't share space] with y.

---

# Intersects

`st_intersects()` is the .hi-medgrn[most general] of confirmation functions
  * returns `TRUE` if x and y geometry share .hi-medgrn[any space], `FALSE` if .hi-medgrn[absolutely no shared space]

Let's use it to retrieve the counties adjacent to Ingham.

```r
ingh_int <- st_intersects(counties, ingham)
class(ingh_int)
```

```
## [1] "sgbp" "list"
```

---

# Intersects

`st_intersects()` is the .hi-medgrn[most general] of confirmation functions
  * returns `TRUE` if x and y geometry share .hi-medgrn[any space], `FALSE` if .hi-medgrn[absolutely no shared space]

By default, `st_intersects()` returns a .hi-medgrn["sparse geometric binary predicate"] (`sgbp`) list
  * 1 list element per geometry in x (wellheads)
  * Each =1 if fully within y (ingham county)
  * =0 if not fully within

---
# Sparse vs. Non-Sparse

Setting `sparse = FALSE` yields the .hi-blue[entire TRUE/FALSE matrix]
* 1 row per row of X
* 1 column per row of Y
* element `i,j` whether geometry row `i` of X intersects geometry row `j` of Y

```r
ingh_int_mat <- st_intersects(counties, ingham, sparse = FALSE)
class(ingh_int_mat)
```

```
## [1] "matrix" "array"
```

```r
head(ingh_int_mat)
```

```
##       [,1]
## [1,] FALSE
## [2,] FALSE
## [3,] FALSE
## [4,] FALSE
## [5,] FALSE
## [6,] FALSE
```

---
# st_filter
We could use index positions and `which` to extract matching geometries from the spatial dataframe...

But `st_filter()` makes the process of retrieving the matches even easier

* returns a .hi-medgrn[spatial dataframe] with the geometries that return TRUE using the chosen predicate
* default predicate: `st_intersects()`
* change predicate with the `.predicate` optional argument

```r
ingh_int_sf <- st_filter(counties, ingham)
ingh_int_sf$NAME
```

```
## [1] "Clinton"    "Eaton"      "Ingham"     "Jackson"    "Livingston"
## [6] "Shiawassee"
```
---
# st_filter

```r
ggplot() +
  geom_sf(data = ingh_int_sf) + 
  maptheme
```

---
# Intersect vs. Within
.hi-purple[Within] is a .hi-purple[stricter] operation than .hi-medgrn[intersection]

* `st_intersects()` returns both the "within" and "overlap" cases

* `st_within()` returns only the "completely within" cases

Let's find which protected areas lie .hi-purple[fully within] Ingham county:

```r
# Verify that CRS is the same as the rail line
stopifnot(st_crs(ingham) == st_crs(whpa))

whpa_within <- st_filter(whpa, ingham, .predicate = st_within)
```

---
# Intersect vs. Within

Comparing .hi-purple[within] to .hi-medgrn[intersects] highlights the difference:

.font90[

```r
whpa_int <- st_filter(whpa, ingham)
ggplot() +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = whpa_int, fill = "green") +
  geom_sf(data = whpa_within, fill = "white") + 
  maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/int_within_plot-1.png" width="60%" style="display: block; margin: auto;" />
]
---

# Within vs. Contains

`st_contains()` is the .hi-turquoise[inverse] of `st_within()`

* `st_contains(x,y)` returns the .hi-turquoise[transpose] of `st_within(y,x)`

```r
whpa_contain <- st_contains(ingham, whpa, sparse = FALSE) # 1 x 2837
whpa_wthn <- st_within(whpa, ingham, sparse = FALSE) # 2837 x 1
identical(whpa_contain, t(whpa_wthn))
```

```
## [1] TRUE
```

---

# Overlaps

`st_overlaps()` yields the protected areas .hi-pink[partially in] Ingham County

```r
ggplot() +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = st_filter(whpa, ingham, .predicate = st_overlaps), fill = "green")  + 
  maptheme
```

---

# Covered by
`st_covered_by` yields the protected areas .hi-slate[completely covered by] Ingham County
* `st_covers` is the inverse

```r
ggplot() +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = st_filter(whpa, ingham, .predicate = st_covered_by), fill = "purple")  + 
  maptheme
```

---

# Touches
`st_touches` yields the protected areas that .hi-green[only touch the boundary] of Ingham County
* Spoiler: there aren't any

.font80[

```r
ggplot() +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = st_filter(whpa, ingham, .predicate = st_touches), fill = "red")  + 
  maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/touches-1.png" width="63%" style="display: block; margin: auto;" />
]
---

# Disjoint
`st_disjoint` returns the elements of x that .hi-grey[don't share space] with y

```r
ggplot() +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = st_filter(whpa, ingham, .predicate = st_disjoint), fill = "green")  + 
  maptheme
```

---
# Spatial Queries: Relationship

We can also **create new geometries** based on the relationship between two spatial objects

* `st_intersection()` creates geometry equal to the .hi-blue[shared portion] of x and y

* `st_union()` combines geometries of x and y into a .hi-green[single geometry]

* `st_difference()` creates geometry equal to the .hi-pink[parts of x] that .hi-pink[don't intersect] y

* `st_sym_difference()` creates geometry equal to the .hi-purple[parts of x and y] that .hi-purple[don't intersect]

* `st_crop()` creates geometry from portion of x that .hi-medgrn[intersects a specified rectangle]

---
# Intersection

`st_intersection()` trims the protected areas to .hi-blue[only the parts contained within] Ingham county

```r
ggplot() +
  geom_sf(data = st_filter(whpa, ingham, .predicate = st_overlaps), fill = "yellow", alpha =0.5) +
  geom_sf(data = ingham, alpha = 0.5) +
  geom_sf(data = st_intersection(whpa, ingham), fill = "blue")  + 
  maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/intersection-1.png" width="50%" style="display: block; margin: auto;" />
---
# Intersection

Note that `st_intersection()` returns attributes from .hi-blue[both intersecting features].
* Mindlessly, without checking whether they still makes sense for the .hi-blue[intersection] geometry.
* Now we have two different sets of area/length variables, plus acre/mile values for **all of Ingham County**

```r
colnames(st_intersection(whpa, ingham))
```

```
##  [1] "OBJECTID"     "WSSN"         "SYSTEM"       "TYPE"         "APPROVAL_D"  
##  [6] "CreatedUse"   "CreatedDat"   "LastEdited"   "LastEdit_1"   "EgleSdeEGL"  
## [11] "ShapeSTAre"   "ShapeSTLen"   "OBJECTID.1"   "FIPSCODE"     "FIPSNUM"     
## [16] "NAME"         "LABEL"        "TYPE.1"       "CNTY_CODE"    "SQKM"        
## [21] "SQMILES"      "ACRES"        "VER"          "LAYOUT"       "PENINSULA"   
## [26] "ShapeSTAre.1" "ShapeSTLen.1" "rail_dist"    "geometry"
```

**Be careful when interpreting attributes after geometry operations.**

---
# Union

`st_union()` combines .hi-green[each geometry] in x with .hi-green[each geometry of y]

* if a protected area falls fully within y, returns y's geometry
* if it falls partially outside, returns the conjoined geometry

.font80[

```r
# calculate union
whpa_union <- st_union(whpa_int, ingham)
# plotting the new, unioned geometry in the fourth row
ggplot() +
  geom_sf(data = whpa_union[4,], alpha = 0.5, fill = "forestgreen")  + maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/union-1.png" width="45%" style="display: block; margin: auto;" />
]

---
# Difference

`st_difference()` returns the .hi-pink[parts] of x (protected areas) .hi-pink[not in]  y (Ingham)

.font90[

```r
ggplot() +
  geom_sf(data = st_difference(whpa_int, ingham), alpha = 0.5, fill = "orangered")  + 
  maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/diff-1.png" width="70%" style="display: block; margin: auto;" />
]

---
# Symmetric Difference
`st_sym_difference()` returns a geometry for each row of x that also includes the .hi-purple[parts of y] (Ingham) that are .hi-purple[not in] x (note the little missing overlapping piece)

.font90[

```r
ggplot() +
  geom_sf(data = st_sym_difference(whpa_int, ingham)[4,], alpha = 0.5, fill = "violet")  + 
  maptheme
```

<img src="data:image/png;base64,#Spatial-Vector_files/figure-html/symdiff-1.png" width="65%" style="display: block; margin: auto;" />
]
---
# Crop
`st_crop()` trims x to a .hi-medgrn[specified rectangle]

* .hi-medgrn[Easiest solution:] use `st_bbox()` as the y object

* .hi-medgrn[Less easy:] use `ymin`, `ymax`, `xmin`, `xmax` arguments to make the box manually

```r
st_crop(whpa_int, st_bbox(ingham))
```

---
# Crop

The cropping is a little hard to see with Ingham county's boundary

Let's convert the bounding box to an sf object so we can plot it instead.

```r
# Create a dataframe with the desired point locations
box <- data.frame(lon = c(st_bbox(ingham)["xmin"], st_bbox(ingham)["xmax"], 
                          st_bbox(ingham)["xmax"], st_bbox(ingham)["xmin"]
                          ,st_bbox(ingham)["xmin"] ), 
                  lat = c(st_bbox(ingham)["ymin"], st_bbox(ingham)["ymin"],
                          st_bbox(ingham)["ymax"],st_bbox(ingham)["ymax"], 
                          st_bbox(ingham)["ymin"])) %>%
  # convert to sf object
  st_as_sf(coords = c("lon", "lat"), crs = st_crs(ingham))%>% 
  # combine the geometry and cast into a single polygon
  summarise((geometry = st_combine(geometry))) %>% 
  st_cast("POLYGON")
```

---
# Crop
  
Plotting the .hi-medgrn[uncropped] protected areas and the bounding box:

.font90[

```r
ggplot() +
  geom_sf(data = box, alpha = 0.5) +
  geom_sf(data = whpa_int, alpha = 0.5, fill = "deepskyblue3") + 
  maptheme
```

]
---
# Crop
  
Plotting the .hi-blue[cropped] protected areas and the bounding box:

.font90[

```r
ggplot() +
  geom_sf(data = box, alpha = 0.5) +
  geom_sf(data = st_crop(whpa_int, st_bbox(ingham)), alpha = 0.5, fill = "aquamarine3")  + 
  maptheme
```

]
---
# Table of Contents:

1. [Spatial Queries: Measurement](#queries_m)

1. [Spatial Queries: Relationship](#queries_r)

1. [Geometric Operations](#geom)

1. [Spatial Joins](#join)